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ABSTRACT 


This paper presents a procedure for modeling structures containing piezoelectric actuators 
using MSC/NASTRAN and MATLAB. The paper describes the utility and functionality 
of one set of validated modeling tools. The tools described herein use MSC/NASTRAN 
to model the structure with piezoelectric actuators and a thermally induced strain to 
model straining of the actuators due to an applied voltage field. MATLAB scripts are 
used to assemble the dynamic equations and to generate frequency response functions. 
The application of these tools is discussed using a cantilever aluminum beam with a 
surface mounted piezoelectric actuator as a sample problem. Software in the form of 
MSC/NASTRAN DMAP input commands, MATLAB scripts, and a step-by-step 
procedure to solve the example problem are provided. Analysis results are generated in 
terms of frequency response functions from deflection and strain data as a function of 
input voltage to the actuator. 


INTRODUCTION 


NASA Langley Research Center, Industry, and Academia have been actively studying 
and developing induced strain actuation devices for aircraft and aerospace applications 
since the late IPSO’s 1 ' 2 . Strain actuation produced as a result of a thermal load and 
piezoelectricity 3 involves straining of components due to temperature changes and 
voltage changes. Piezoelectric materials such as Lead Zirconate Titanate (PZT) ceramics 
when subjected to an electric field produce mechanical strain or alternately generate an 
electric charge when subjected to a mechanical strain. This property gives piezoelectric 
materials the ability to act as actuators or sensors. Using piezoelectric actuators and 
sensors to form self-controlling and self-monitoring systems to improve performance of 
aircraft and space structures has attracted interest in the research community. 


Numerous studies have been completed on analyses and models for piezo-electrically 
controlled structures. Some of these studies include: a high-order theory to model 
composite laminates with surface bonded or embedded piezoelectric sensors by Seely 4 ; a 
three-dimensional finite element code to analyze the mechanical-electrical response of 
laminated composites containing distributed piezoelectric ceramics developed by Sung 
Kyu Ha 5 ; the use of classical laminate theory to estimate the through-the-thickness strain 
distribution of composite laminates with embedded actuators by Crawley 6 and others 7 ' 8 . 
Although those analytical techniques showed good correlation with experimental data, 
few have been incorporated into commercially available codes for general use. 


Due to the increased interest in the design of complex structures with piezoelectric 
actuators and the need for fast and simple implementation of piezoelectric control 
systems, technology developers are beginning to incorporate or provide the tools to 
create piezoelectric elements. Freed 9 developed one and two-dimensional finite 
elements, which include piezoelectric coupling, for integration into MSC/NASTRAN, 
however, this option is not currently available. Hauch 10 investigated using ABAQUS 
electromechanical-coupled finite elements and superelement capabilities for modeling 
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piezoelectric actuators have allowed a number of viable analytical and numerical tools. 
Reaves 11 developed benchmark test article structures to validate techniques for modeling 
structures containing piezoelectric actuators using commercially available Finite Element 
Analysis (FEA) packages. This document presents detail information on the 
implementation of one approach described in reference 11. The use of these tools is 
presented using a cantilever aluminum beam with surface mounted piezoelectric actuators 
as a sample case. 

The method described herein use MSC/NASTRAN to model the structure with 
piezoelectric actuators and a thermally induced strain to model expansion/contraction of 
the actuators due to an applied voltage field. To reduce the number of structural modes 
needed for an accurate solution, Ritz vectors 12 are appended to the structural modes. 
MATLAB scripts are used to assemble the dynamic equations and to generate frequency 
response function at the locations of interest. Software in the form of MSC/NASTRAN 
DMAP input commands, MATLAB scripts, and a step-by-step procedure to solve the 
case history example are provided. Frequency response functions from deflection and 
strain responses to applied input voltage are used for evaluation of the methodology. 


PROBLEM FORMULATION 


Modeling of a structure with piezoelectric actuators using finite elements can be 
performed like any conventional finite element model formulation. However, to include 
piezo effects, elements that expand or contract as voltage is applied, the thermal strain 
capability within the element formulation needs to be used to evaluate voltage induced 
internal forces. The governing equation of motion is written as 

Mx + Kx = f p (1) 

where f represents internal loads due to actuation, M is the mass matrix, K is the 

stiffness matrix, and jc contains node displacements. Induced forces are computed by first 
defining an artificial coefficient of thermal expansion a =d 31 /t using the piezoelectric 

strain to voltage parameter d 3l and the piezo material thickness t. The second step is to 
apply a unit change in temperature (or voltage) in NASTRAN to evaluate the resulting 
deformation T r . This deformation vector T r is referred to as a Ritz vector, one Ritz vector 
is required for each independently actuated piezo. Internal loads are now related to 
displacements and voltages by / = KT r Av . In the reduced order model formulation, 

mass normalized eigenvectors <f> from an eigenvalue solution of Eq. (1) and a set Ritz 
vectors T r can be used as a set of basis vectors for model reduction. With these basis 
vectors a new coordinate transformation can be defined 

*■1® r ']{zJ (2) 

to transform Eq. (1) to 


2 



Av 


(3) 


' 1 

M r - 

f z } 

< ^ + 

A 


H- 



M rr_ 

L^rJ 

X 


UJ 

Krr. 


where A is a diagonal matrix containing the eigenvalues to the homogeneous form of Eq. 
( 1 ), I is an identity matrix, and all the other matrices are partitioned consistent with the 
number of basis vectors used. For example, using m eigenvalues and r piezoelectric 
actuators, the transformation matrix in Eq. (2) has m+r columns and n rows, where n is 
the number of degrees-of-ffeedom in the model. Similarly the transformed displacement 
vector z has m and z r has r rows. Often, Eq. (3) is sufficient to conduct simulation 

studies, however, in this form the equation corresponds to a stiff system of differential 
equations because the transformed mass matrix is ill conditioned. Ill conditioning in this 
particular case is due to negligible mass values associated with the Ritz basis vectors. To 
overcome this problem, one can decompose the mass matrix using singular value 
decomposition such that 







(4) 


where U contains the left singular vectors, V contains the right singular vectors, and cr 
are nonzero singular values. In the numerical implementation of the procedure it is not 
likely that some of the singular values will be identically equal to zero but values below a 
certain threshold are considered zero. Furthermore, the number of singular values with 
numerical values near zero equals the number of piezoelectric actuators placed in the 
model. Using the right singular vectors, a new transformation matrix can be defined as 


z 



(5) 


Substituting Eq. (5) into Eq. (3) yields a transformed set of equations 
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where the bottom set of equations correspond to an algebraic set of equations. Separating 
the algebraic equation from the differential equations, the resulting properly conditioned 
set of matrix differential equation is 

crrj + (* - k r ky r )n = (b- kX %) ^ (7) 


Time integration of Eq. (7) provides values for 77 and the solution for 77 ,. is 

ri r = ~k^k T r 77 + k^b r Av ( 8 ) 
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To recover the response in physical coordinates, results from Eqs. (7) and (8) are 
substituted into Eqs. (5) and (2). For MATLAB implementation all the equations are 
transformed into state space form and then integrated using the MATLAB built-in 
integration software. 


EXAMPLE PROBLEM 


A cantilevered aluminum beam 2.75” x 16” x .04”, Figure 1, with one piezoelectric 
actuator bonded near the root is used in this study. A Flex-Patch piezoelectric actuator, 
developed and fabricated at NASA LaRC, is selected for the application. The Flex-Patch 
consists of a 3”xl.75”x.008” Morgan Matrox PZT-5 A 13 piezoceramic encapsulated using 
a polymer film. The PZT-5A piezoceramic mechanical and electro-mechanical 
properties are listed in Table 1. MSC/NASTRAN is used to model the structure; Figure 2 
shows a finite element representation of the cantilevered beam indicating the location of 
the piezoelectric actuator and strain sensor output used in the present study. The model 
contains 396 quadrilateral elements and 450 nodes. PCOMP cards in MSC/NASTRAN 
are used to specify the properties of the composite lay-up for the actuator and aluminum 
substrate. 


ANALYTICAL TOOLS IMPLEMENTATION 


The procedure to generate the analytical model of the cantilever beam involves the 
following steps: 1) construction of a FEM based on the test article and Ritz vectors 
generation, 2) eigenvectors generation and mass and stiffness matrices calculation to 
assemble the reduced order model, 3) execution of MATLAB scripts to assemble the 
dynamic equation and to generate frequency response function at the locations of interest. 
These steps are described in more detail in the following sections. 

Structural FEM and static Ritz vector generation - MSC/NASTRAN offers no 
piezoelectric coupled-field elements capability from which to model smart structures 
directly. Rather, the analogy between piezoelectric strain and thermally induced strain, is 
used to allow temperature changes to model piezoelectric voltage actuation. Piezoelectric 
coefficients characterizing the actuator are input as thermal expansion coefficients 
(CTE’s) associated with standard elements. In this study the model treats both the 
actuator and the structures substrates as plies of an integrated laminated plate. PCOMP 
cards in MSC/NASTRAN are used to specify the properties of the composite lay-up and 
the applied voltage is modeled as a thermal load. TEMP cards identify the locations at 
which thermal loads (voltage) is applied. MAT1 cards representing the actuators include 
non-zero thermal expansion coefficients. The thermal expansion coefficients in the 
remaining structure must be zero to ensure thermal expansion loads are generated only 
at the piezoelectric actuator locations. Appendix 1 shows a listing of the 
MSC/NASTRAN input data deck ‘almq2.dat’, for the cantilevered beam. For brevity, 
MSC/NASTRAN grid points and element information are not included in the listing. The 
execution of this file will generate the Ritz vector (static deflection and/or strain data for 
the beam) using a thermal load equivalent to 1 volt. The DMAP sequence ‘static.dmap’ 
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incorporated in the MSC/NASTRAN input data deck stores the static solution vector in a 
binary formatted output file ‘S101SV1.OUT4'. In applications with multiple actuators 
each actuator load and output request are to be included in individual SUBCASES. 

Eigenvectors and reduced model generation - After computing Ritz vectors, a general 
eigenvalue/eigenvector MSC/NASTRAN solution containing special structural matrix 
transformation routines (DMAP) combine the Ritz vectors with eigenvectors and 
calculates the mass and stiffness matrices in Eq. (3) needed to assemble a reduced order 
model. NASTRAN commands for the cantilever beam example are listed in the file 
named ‘almq2e.dat’ and the DMAP sequence is listed in ‘genkm^Z.dmap’ (Appendix 
2). The ‘almq2e.dat’ file uses the file ‘S101SV1.OUT4’ as input and generates two 
ASCII formatted output files: ‘almq2e.asc’ and ‘almq2e.pch’ containing the reduced 
mass/stiffness matrix and the displacements and/or strains for the structural modes 
including the Ritz vectors, respectively. Output observation points are for displacements 
at nodes 348 and 349 and strain for element 116. To facilitate importing the matrices to a 
MATLAB environment the file ‘almq2e.asc’ needs to be edited to remove lines starting 
with the one containing “eigensum” up to but not including the line containing “genk”. 
The first line in the file should contain ‘genk’ and must not have blank lines in the 
beginning of the file. 

Dynamic equation assembly and results generation - MATLAB scripts are required to 
assemble the dynamic equations and generate frequency response functions at the 
observation points. The MATLAB script ‘makemodel.m , used in the cantilever beam 
case study is listed in Appendix 3. The script creates space state model matrices based on 
the number of inputs (number of actuators) and MSC/NASTRAN output request (element 
strain and/or nodal displacement). The script uses as input the files generated and edited 
in the previous step: ‘almq2e.asc’ and ‘almq2e.pch’. The MATLAB scripts are problem 
specific depending on the type of output requests, but are easy to modify. 

A flow chart summarizing the complete procedure to solve the cantilevered aluminum 
beam problem is shown in Figure 3. 


EXAMPLE PROBLEM RESULTS 


The frequency response functions for input voltage, tip displacement, and strain gage data 
are generated from analysis. NASTRAN analysis results using the modeling technique 
described herein are shown in Figure 4 for the first four beam/actuator mode shape 
deformations with corresponding frequencies of 5.68, 33.59, 60.24 and 91.21 Hz. To 
examine the input/output relationship of the system with the actuator. Figure 5 shows a 
comparison of the frequency response function magnitude (top) and phase (bottom) from 
the beam tip displacement to the piezo-actuator input voltage. Furthermore, Figure 6 
shows results from the nearly collocated strain gage. 
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SUMMARY 


A procedure for modeling structures containing piezoelectric actuators using 
MSC/NASTRAN and MATLAB is presented. The procedure describing how to use 
these tools is presented using a cantilever aluminum beam with surface mounted 
piezoelectric actuators as example problem. Software in the form of MSC/NASTRAN 
DMAP input commands, MATLAB scripts, and a step-by-step procedure to solve the 
example problem are provided. Analysis results are generated in terms of frequency 
response functions from deflection and strain data to input voltage to the actuator. 
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Table 1. Material properties of PZT-5A piezoceramic and T300/976 graphite/epoxy 

composite 



PZT-5A 

T300/976 

Modulus of elasticity 
(lbf/in 2 ) 



Ei 

1.0E+7 

2.17E+7 

e 2 


1.305E+6 

Poisson’s ratio 



V 

0.3 

0.3 

Shear Modulus (lbf/in 2) 



Gi2 

3.82E+6 

1.03E+6 

Giz 


1.03E+6 

G 2 z 


363600 

Density, (lbf- sec 2 /in 4) 



P 

7.16E-4 

1.49E-4 

Piezoelectric constant, 
(in/Volt) 



cbi 

6.73E-9 


Electrical permittivity, 
(farads/in) 



s! 1 

4.2E-10 
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Actuator 



Figure 1. Cantilevered aluminum beam with piezoelectric actuator. 



El 16 - strain sensor output location 
Nodes 348, 349 - displacement output location 


Figure 2. Aluminum beam with piezoelectric actuator finite element mesh. 
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Run MATLAB script makemodel.m 
Creates space state model 


Figure 3. Steps to generate analytical model of aluminum beam with piezoelectric 

actuator. 
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Figure 4. Mode Shapes of aluminum beam with piezoelectric actuator. 



Phase 

(Degrees) 



Frequency (Hz) 


Figure 5. FRF of aluminum beam tip displacement as a function of piezoelectric actuator 

input voltage. 
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Appendix 1 


The following NASTRAN input data deck listing is contained in the almq2.dat file. The 
deck is set up to generate the Ritz vector corresponding to a single piezoactuator mounted 
on the structure. TEMP cards identify the locations at which thermal loads (voltage) is 
applied. The DMAP static.dmap incorporated in the data deck by the include" 
command stores static solution vectors in the output file ’S101SV1.OUT4’ specified in 
the ASSIGN statement. 

almq2.dat data deck: 

ASSIGN OUTPUT4 =, S101SVl.OUT4' UNIT=15 DELETE 
INIT MASTER(S) 

ID PIEZOELECTRIC QUAD ELEMENTS 
SOL 101 
$ 

compile sestatic,souin=mscsou,nolist,noref$ 
include 'static.dmap' 

$ 

TIME = 60 
CEND 

TITLE = NASTRAN FILE TRANSLATOR - UNITS = SI 
$ GLOBAL CASE 
SPC = 1 

TEMP = 1 

$ Define set of elements for strain output 

SET 4 =116 

$ Define set of nodes for displacement output 

SET 5 =348,349 
DISPLACEMENT(PRINT)=5 
STRAIN(PRINT,FIBER)=4 
S 

BEGIN BULK 
$ 

PARAM AUTOSPC YES 
PARAM POST -2 
PARAM,K6ROT, 1 00. 

GRID 1 0 7.14E-3 5.75E-2 0.00000 

GRID 2 0 7.14E-3 4.87E-2 0.00000 

GRID 3 0 7. 1 4E-3 3.98E-2 0.00000 


CQUAD4 396 2 552 513 519 


0 

0 

0 


555-90.0000 
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$ Property card definition for aluminum beam substrate 

PCOMP 2-5.08E-4 0.00000 0.00000 0.00000 + 

+ 2 1.27E-5 0.00000 YES 2 9.91E-4 0.00000YES + 

+ 2 1.27E-5 0.00000 YES 


$ Property card definition for aluminum beam substrate and piezoactuator lay-up 


PCOMP 

5-5.08E-4 0.00000 

0.00000 0.00000 

+ 

+ 

2 1.27E-5 0.00000YES 

2 9.91E-4 0.00000 YES 

+ 

+ 

2 1.27E-5 0.00000 YES 

3 1.27E-5 0.00000YES 

+ 

+ 

3 2.54E-5 0.00000YES 

4 2.03E-4 0.00000YES 

+ 

+ 

3 2.54E-5 0.00000 YES 

3 1.27E-5 0.00000 YES 



$ Aluminum material properties 

MAT1 2 7.2E+10 0.33333 2758. 0.0 


SKapton film material properties 

MAT1 3 2.76E+9 0.45000 1358. 0.0 


$Piezo ceramic material properties, non-zero CTE 

MAT1 4 6.9E+10 0.31000 7700. 8.42E-7 

$ 

SPC1 1 123456 504THRU 506 

SPC1 1 123456 519THRU 523 

SPC1 1 123456 554THRU 555 

$ Unit thermal load (voltage) applied to actuator 


TEMP 

1 

1 

1.000 

2 

1.000 

3 1 

1.000 

TEMP 

1 

4 

1.000 

5 

1.000 

6 1 

1.000 

TEMP 

1 

7 

1.000 

8 

1.000 

9 i 

1.000 

TEMP 

1 

10 

1.000 

11 

1.000 

12 

1.000 

TEMP 

1 

13 

1.000 

14 

1.000 

15 

1.000 

TEMP 

1 

16 

1.000 

17 

1.000 

18 

1.000 

TEMP 

1 

19 

1.000 

20 

1.000 

21 

1.000 

TEMP 

1 

22 

1.000 

23 

1.000 

24 

1.000 

TEMP 

1 

25 

1.000 

26 

1.000 

27 

1.000 

TEMP 

1 

28 

1.000 

29 

1.000 

30 

1.000 

TEMP 

1 

31 

1.000 

32 

1.000 

33 

1.000 

TEMP 

1 

34 

1.000 

35 

1.000 

36 

1.000 

TEMP 

1 

37 

1.000 

38 

1.000 

39 

1.000 

TEMP 

1 

40 

1.000 

41 

1.000 

42 

1.000 

TEMP 

1 

43 

1.000 

44 

1.000 

45 

1.000 

TEMP 

1 

46 

1.000 

47 

1.000 

48 

1.000 

TEMP 

1 

49 

1.000 

50 

1.000 

51 

1.000 

TEMP 

TEMPD 

ENDDATA 

1 52 

1 0.0 

1.000 

53 

1.000 

54 

1.000 
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Below is a listing of the ‘static.dmap’ DMAP: 

$compile sestatic, souin=mscsou, nolist, noref $ 

$ 

$ Following DMAP was coded for MSC/NASTRAN and verified to 
$ work in version 2001 .0.1 for the current application 
$ MSC says this will work for versions 68.2,69, 69.1 or 70 
$ This DMAP stores static solution vectors in an output file 
$ 'S101SV1.OUT4’ specified in the ASSIGN statement. 

$ If the FEM does not contain rigid body modes (or mechanisms), 

$ then the L-set d.o.f. are written to the output file. 

$ If there are rigid body modes (a SUPORT card should be included 
$ in the Bulk Data deck), then the T-set d.o.f. are written 
$ to the output file. 

$ 

alter 'CALL SUPER3’ (1,0) 

$alter 685 $ after static solution SSG3 
$ If there is no R-set, write the L-set vectors to file 
if (norset < 0) then $ 
output4 ul„„//-1/15/1 $ 

$ Print L-set solution vector 
$ matprn ul II $ 
endif $ 

$ 

$ If the R-set exists, write the T-set vectors to file 
Salter 822 $ After SDR1 module 
if (norset <>-1) then $ 

$ Partition the G-set solution, UGVS, into the T-set vectors 
if (nogset = notset) then $ no partitions if g-set=t-set 
output4 ug„„//-1/15/1 $ 

$ matprn ug II $ 
else $ 

upartn uset,ug/tugvs„,/'g7't'/'o71 $ 

$ matprn tugvs// $ T-set solution displacement vector 
output4 tugvs„„//-1/15/1 $ 
endif $ end of "if (nogset=notset)" 
endif $ end of "if (norset <> -1 )" 

$ 

endalter 
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Appendix 2 


The following NASTRAN input data deck listing is contained in the almq2e.dat file. The 
executive deck is set up to generate vibration modes. The DMAP sequence 
genkm_v2.dmap incorporated in the data deck by the ‘include’ command combines the 
Rjtz vectors and eigenvectors and calculates the mass and stiffness matrices needed to 
assemble a reduced order model. The almq2e.dat file uses the file S101SV1.OUT4 
generated by executing almq2.dat as input and generates two ASCII format output files: 
almq2e.asc and almq2e.pch containing the reduced mass and stiffness matrix and the 
displacements and/or strains for the structural modes including the Ritz vectors 
respectively. 

ASSIGN INPUTT4— S101SV1.OUT4' UNIT=15 

ASSIGN OUTPUT4-almq2e.asc' UNIT=16 FORM=FORMATTED DELETE 
ID ALBEAM WITH ACTUATOR,PROB 
SOL 103 

include 'genkm_v2.dmap’ 

TIME = 60 
CEND 

TITLE = NASTRAN FILE TRANSLATOR - UNITS = SI 
$ GLOBAL CASE 
SPC = 1 

METHOD=25 

$ Output requests are for displacements at nodes 348 and 349 
$ and strain for element 1 16. 

S 

SET 4 =116 
SET 5 =348,349 
DISPLACEMENT (PUN C H)=5 
STRAIN (PUN CH,FIBER)=4 

$ 

PARTN=5 

$ 

BEGIN BULK 
EIGRL,25,„10 
$ 

PARAM,K6ROT, 1 00. 

PARAM AUTOSPC YES 
PARAM POST -2 

GRID 1 0 7. 14E-3 5.75E-2 0.00000 0 


The rest of the data deck is identical to the listing included in ‘almq2.dat’ in Appendix 1. 
TEMP cards are ignored in an eigenvalue analysis. 
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The following is a listing of the ‘genkm_v2.dmap’ DMAP: 

$MSC/NASTRAN DMAP 

$ this dmap converted to V69 rjb 1/14/98 

$This dmap was verified to work for version 2001 .0.1 

WYTt t TTT YtTtY TTTTTTYttTTTTTt YTTt t T t t TTTTTT 

$ This DMAP accepts user selected static solution vectors from an 
$ input file. The projection of these vectors onto a subspace 
$ perpendicular to the subspace of mass weighted normal modes 
$ is computed. These new static vectors are appended to the 
$ calculated normal modes, and generalized stiffness and mass 
$ matrices are formed. The static partition of these generalized 
$ matrices are placed on an ASCII output file. Also written to 
$ the ASCII output file are: a normal modes summary matrix, 

$ and user selected static vectors in external sort order. 

$ 

malter 'CALL.*IFPL.*DMINDX' $ after dmindx is formed 
call dbstore dmi, dmindx, „//99/99/' 70 $ 

$ 

compile moders souin=mscsou, nolist, noref $ 

alter 'LAMA,OEIGS'{1 ,1) $ after REIGL and OFP modules 

$ NEIGV = number of eigenvalues computed 

$ NSV = number of columns in static vector matrix STATV 

$ STATV is either an L-set or T-set dimensioned static vector matrix 

inputt4 /statv„„/1/15/-1/1 $ read matrix of static vectors 

$ 

$ Matrix SSV is an optional input. Matrix SSV is used to select 
$ which static vectors will be used in the following analysis. 

$ If matrix SSV is not input, then all the static vectors in the 
$ inputt4 file will be used in the analysis. 

$ 

$ DMIIN module is required to read the DMI bulk data entries. 

$ The name of the input DMI data block is SSV. 

$ 

type parm„i,n,indxok 

call dbfetch /dmi, dmindx,,, /99/99/0/0/s,indxok 
message II' indxok = Vindxok $ 

dmiin dmi,dmindx/ssv, ,„/s,n,yesv $ 

if (yesv) then $ If matrix SSV is input perform following operations 
$ Postmultiply STATV by SSV to select the desired static vectors. 

$ The resulting vectors are stored in DSV. 
mpyad statv,ssv,/dsv $ 

$ matprn dsv// $ 

$ Let DSV overwrite STATV 

equivx dsv/statv/-1 $ DSV is renamed to STATV 

$ matprn statv// $ 

endif $ Above "IF (yesv)" group only used if SSV exists in bulk data 
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$ 

paraml statv//'trailer'/1/s,n,nsv $ nsv=number of columns in STATV 
$ message //'number of static vectors=7nsv/' modes='/neigv $ 

$ matprn statv,phix,kxx,mxx// $ 

param //'add7s,n,nritz/neigv/nsv $nritz=number of modal+static vectors 
$ Form the following summation - 
$ [sproj]=[statvj - [phix]*[phix]'*[mxx]*[statv] 

Ssmpyad phix.phix , mxx.statv, ,statv/sproj/4/- 1 /I /// 1 $ 

$ matprn sproj II $ 

$ Append columns of STATV to columns of PHIX to form PHINEW 
matgen ,/cpvec/6/nritz/neigv/nsv $ column partitioning vector 
merge phix„statv„cpvec,/phinew/1 $ 

$ matprn phinew II $ 

$ 

$ Compute system generalized stiffness and mass matrices 
smpyad phinew, kxx,phinew,„/genk/3/1///1////6 $Symmetric stiffness matrix 
smpyad phinew.mxx, phinew,, ,/genm/3/1///1////6 $ Symmetric mass matrix 
$ matprn genk,genm// $ 

$ Extract the (NSV by NSV) static vector partition from the generalized 
$ stiffness and mass matrices 
Spartn genk,cpvec,/,„sgenk/-1 $ 

$partn genm,cpvec,/,„sgenm/-1 $ 

$ matprn sgenk,sgenm// $ 

lamx, ,lama/eigsum/-1 $ Convert LAMA from table to matrix format 
$ Expand LAMA to include the static vectors 
matgen ,/pvfm/6/5/2/2/1 $ 
partn eigsum,pvfm,/„fandm,/1 $ 

$ matprn fandm// $ 
matgen ,/xcol/6/3/1/1/1 $ 
matgen ,/xrow/6/nritz/neigv/nsv $ 

merge fandm„„xcol,xrow/mlama/1 $ expanded transposed LAMA matrix 
$ matprn mlama II $ 
tmsp mlama/tlama $ 

lamx tlama./lamnew $ Convert LAMA matrix back into table form 
equivx lamnew/lama/-1 $ Store new table in LAMA 
equivx phinew/phix/-1 $ Store augmented eigenvectors in PHIX 
neigv = nritz $ NEIGV=no. of modes + no. of static vectors 
$ 

output4 eigsum,genk,genm„//-1/16/0 $ 

$ 

V'P V 4)4)vPvPvbvP vbyD wOwCDuJw wu)u) J) \J)vPvPU) 4) \J)\P4)\J) 

compile sedrcvr nolist noref 
$alter 1157,1157 $ Before SDR2 
alter 'SDR2'(1,-1) $ 

$ Convert displacements to the basic coordinate system 
vecplot ug,bgpdts,eqexins,cstms,casedr,„/ugvb/0/0/1////// $ 

$ 
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$ Convert from internal to external sort 
matgen eqexins/intext/9//lusets $ 
mpyad intext, ugvb,/xgall/1 $ 

$matprn xgall// $ 

$ 

$ Using Case Control input PARTN, create a partitioning vector for 
$ the grids whose output is desired. The partitioning vector has 
$ G-set dimension, and a value of 1 .0 associated with the d.o.f. 
matmod eqexins,uset,sils,casedr,,/ucpv,/17//1 $ 

$ matprn ucpv// $ 

$ Convert partitioning vector to external sort order 
mpyad intext, ucpv, /xucpv/1 $ 

$ matprn xucpv//$ 

$ Apply partitioning vector to extract desired grids 
partn xgall„xucpv/xcomp,xdgp„/1 $ 

$ matprn xdgp II $ 

$ 

$ Output in ASCII format 

$ Note: user must include in File Management Section the next line 
$ assign output4='ritzv.asc' unit=16form=formatted 
output4 xdgp„„//0/16/0 $ 
endalter $ 

(p fl* (r (p (p tp |p (p IP |P ^ ^ y^ |p IP IP |P y^ y^ y^ y y^ |p (p y^ |P (P |P ^ ^ y^ y* y^ ^ ^ ^ |P |P (P |p IP 

^P *p J> w J ) J) ^p ^p ^p ^p ^p ^p ^p ^p ^p ^p ^p jn yP Cp \P%p \P ^p \P p ^p ^p u) ^p ^p ^P ^p ^p r p ^p ^p ^p u) ^p ^p ^p ^p ^p ^p ^p ^p 
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Appendix 3 


The following MATLAB script is used to create space state model matrices based on the 
number of inputs (number of actuators) and MSC/NASTRAN output request (element 
strain and/or nodal displacement). Caution: The MATLAB scripts are problem specific 
depending on the type of output requests, but easy to modify. 


makemodel.mat 

0 ^******** **************************************** ********************** 

% Script file MAKEMODEL reads NASTRAN Files * 

% with the reduced mass, stiffness, mode shapes, and Ritz vectors * 

% * 

% A state space model of the system is constructed with input in terms * 

% of voltages to the piezoactuators and output as desired * 

% * 

% Created by L.G. Horta and M.C. Reaves * 

% NASA Langley 8-01-01 * 

% * 

% 

Fname2=input(' Enter NASTRAN ascii FileName for reduced mass and stiffness (CR=default):= Vs'); 
if isempty(Fname2); Fname2- almq2e.asc';disp(Fname2);end 


[FID,MESSAGE] = fopen(Fname2,'rt' ); 

[Stiff,Mass]=readMK(FID); % Reading reduced mass and stiffness matrices 

Fname2=input(' Enter NASTRAN Punch FileName (CR=default):= \V); 
if isempty(Fname2); Fname2- almq2e.pch’;disp(Fname2);end 
[FID,MESSAGE]= fopen(Fname2/rf ); 

nq=input(' Number of Output nodes (CR=default):= '); 
if isempty(nq); nq=2;disp([ r nq = ' num2str(nq)]);end 
nm=input(' Number of Modes in File (CR=default):= '); 
if isempty(nm); nm=l l ;disp(f' nm = ' num2str(nm)]);end 
ne=input(’ Number of Strain ELEMENTS (CR=default):= '); 
if isempty(ne); ne=l;disp([' ne = ' num2str(ne)]);end 
Nplies=input(’ Number of Plies per element (CR=default):= '); 
if isempty(Nplies); Nplies=3;disp([ , Nplies = ' num2str(Nplies)]);end 
PlyNum=input(' Ply number where sensor data is measured (CR=default):= '); 
if isempty(PlyNum); PlyNum=3;disp([' PlyNum = ' num2str(PlyNum)]);end 
NR=input(’ Enter Number of Ritz- Vectors (CR=default):= '); 
if isempty(NR); NR=1 ;disp([' NR = ’ num2str(NR)]);end 

[eigen,Shape]=readShape(FID,mn,nq); % read node information 
f Strain]=readply(FLD,Nplies,ne,PlyNum,nm); % read element information 
fclose(FID); 


o^*********************** ******************************** ******* ****** 

% End of Read Punch File 

% 

% construct State Space Model 

% 
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ShapeStrain=[Shape; Strain]; % Stack output matrices 

[AJB„CD,ModelDesc]=ConstructABCD70301(NR,Mass,Stiff,ShapeStrain); 

% 

% Output order is defined by NASTRAN File node order followed by element info 

% 

% For Example: Node =[ lx ly lz 2x 2y 2z ...nqx nqy nqz Element 1 ....] 

SYS=ss(A,B,C,D); 

[Noutjunk]=size(C); 

Nfpts=1600; 

ffeq=linspace( . 1 ,400,Nfpts)’; 

System=ss(A,B,C,D); 

[mag,phase]= : bode(System,freq*2*pi); 

Mag=reshape(mag,[Nout*NR Nfpts]); % Stored Columnwise 

Phase=reshape(phase,[Nout*NR Nfpts] ); 
for j=l:NR 
for l=l:Nout 

subplot(2, f 1 ),semilogy(freq,Mag((j- 1 )*Nout+l,:)'); 
xlabelf Freq (Hzy);ylabel(' Mag’); 

title([’ Input No. := 1 num2str(j) ' Output No. := ’ num2str(l)]) 
subplo^ij^^lo^freq^PhaseCO-l^Nout+l,:)’); 
xlabelf Freq. (Hz)’); ylabelf Phase'); 
dispf Hit Carriage return to see next plot’) 
pause 
end 
end 


readMK.mat 

function [ Stiff Mas s]=readMK( FID) 

% 

% Function to read ASCII File from NASTRAN 
% with reduced MASS and Stiffness Matrices 
% 

% The File must be 8 1 columns wide 
% 

% [Stiff,MASS]=readMK(FID) 

% 

% Input Information 

% FID = NASTRAN ASCII file name pointer (81 col. wide) 
% 

% Output Information 
% Mass = Reduced Mass Matrix 
% Stiff = Reduced Stiffness Matrix 
% 

% Written by Lucas G. Horta 01252001 

if FID >= 0; 

LINE=fgets(FID); 

NN=str2num( L INE( 1 :20)); % Two Elements NCOL, NROW 

Mass=zeros(NN(2),NN( 1 )); 

Stiff=zeros(NN(2),NN( 1 )); 
for j=l:NN(l); 

LINE2=fgets(FID); 
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NN2=str2num(LINE2( 1 :20)); % Two elements COL, RowPos 

Mat=[]; 

TT=fscanf(FID.’%g',NN( 1 )); 

TSpace 1 =fgets(FID); 

Stiff(:j)=TT; 

end 

LINE2=fgets(FID); 

LINE3=fgets(FED); 

LINE4=fgets(FID); 
for j=l:NN(l); 

LINE2=fgets(FID); 

NN2=str2num(LINE2( 1 :20)); % Two elements COL, RowPos 

Mat=[]; 

TT=fscanf( FID,’%g' ,NN ( 1 )); 

TSpace=fgets(FID); 

Mass(: j)=TT; 
end 

fclose(FID); 

else 

disp(' File not found’) 
end 

readShape.mat 

function [eigen,Shape]=readShape(FID,nmodes,nq) 

% 

% Function to read PUNCHFILE from NASTRAN 
% with shape information 
% 

% The File must be 8 1 columns wide 

% 

% [eigen,Shape]=readShape(FID,nmodes,nq) 

% 

% Input Information 

% Filename = NASTRAN Punch file name (ASCII) (81 col. wide) 
% nmodes = Number of modes 
% nq = Number of output nodes 

% 

% Output Information 

% eigen = Eigenvalue squared (Rad/sec) 

% Shape = Shape Vector size (ngpt*6 x nmodes) 

% 

% Written by Horta-Reaves 7-03-01 
if FID >= 0; 
for ijk=l:nmodes 
forjj=l:6 

LENE=fgets(FED); %Header 

end 

LINE7=fgets(FID); % Eigenvalue ED 

eigen(ijk)=str2num(LINE7( 14:30)); 
tem=[]; 
for 1=1 :nq 
LINE 1 =fgets(FID); 

[nn,mm]=size(LINE 1 ); 

nYV7=ctr?rmmU IXJF 1 Ofl*mm ^ V 

“ ^ 
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pxyz=pxyz(l:3); 

LINE3=fgets(FID): 

% Rotation DOF removed 

% [nn,mm]=size( LINE3 ); 

% txyz=str2num(LINE3(20:mm)); 

% txyz=txyz(l:3); 

% tem=[tem pxyz txyz]; 

tem=[tem pxyz]; 
end 

Shape(:,ijk)=tem'; 

end 

else 

dispf File not found’) 
end 


readply.mat 

function [Strainl=readply(FID,Nplies,ne,PlyNum,nm) 

% 

% Supporting file to read NASTRAN punch file info from composite elements 
% 

% [Strain]=readply(FID,Nplies,ne,PLyNum,nm) 

% 

% Input 

% FID = File handle 
% Nplies = Number of plies 
% PlyNum = Ply number where output is recovered 
% ne = Number of strain elements 
% nm Number of modes 
% 

% Output 

% Strain^ [ne x nm] with strain influence in the Normal- 1 Direction for all elements 

% 

% 


% Read strain on PLyNum NORMAL- 1 Direction 
% Created 8-01-01 by L. Horta 
Strain==zeros(ne,nm); 

Numberoflines^fPlyNum- 1 )*4; 
for ijk=l:nm 
forjj=l:8 

LINE=fgets(FID); %Header 

end 

tems=[]; 
for 1=1 :ne 

% Read strain on PlyNum, NORMAL- 1 Direction 
if Numberoflines > 0; 

for zl=l .Numberoflines 
LINE=fgets(FID); 
end 
end 

LINE 1 =fgets(FID); 

[nn,mm]=size(LINEl); 
exy=str2num(LINE 1 (20:mm)); 

% Select NORMAL- 1 Direction 
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% To get normal-2 change the number 2 to 3 (See NASTRAN Output convention) 
exy=exy(2); 
tems= [terns exy]; 

for z2=l:(Nplies-PlyNum)*4+3 

LINE=fgets(FID); % Remaining plies info not needed 

end 
end 

StrainOJjk^ems 1 ; 


ConstructABCD70301.mat 

function [A,B,CJXModelDesc]=ConstructABCD70301(NR,Mass,Stiff,Shape) 

% 

% Function to construct A,B,C,D, matrices for a piezoelectric actuated 
% system 
% 

% [A,B,C,D]=ConstructABCD70301(NR,Mass,Stiff,Shape) 

% 

% Input Information 

% NR = Number of Ritz Vectors 

% Mass = Mass Matrix 

% Stiff = Stiffness Matrix 

% Output Information 

% eigen = Eigenvalue squared (Rad/sec) 

% Shape = Shape Vector size (Number of Outputs x nmodes) 

% ModelDesc = Text Describing output definition 

% 


% Written by Horta-Reaves 7-03-01 (Mod 8-01-01) 


% 

% construct State Space Model with Ritz vector and mode shape info 

% 

[NM,Njunk]=size(Mass); 


% 

% Reduction of Algebraic equation for Ritz Vectors 

% 

I1=[1:NM-NR]’; 

I2=[NM-NR+1 :NM] r ; 

[u,s,v]=svd(Mass); 

Index=abs(di ag( s))> 1 e- 1 0; 

NMnew=sum(Index); 

Knew=u'* Stiff 5 v; 

Bnew=u f * Stiff( :,I2); 

Kll=Knew(IUl); 

K12-Knew(Il,I2); 

K22=Knew( 12,12); 

bl=Bnew(Il,:); 

b2=Bnew(I2,:); 

K22in=inv(K22); 

CP=Shape*v; 

C1_P“CP(;,I1); 
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C2_P=CP(:,I2); 

% Mod 3-28-01 
[Nsl Junk]=size( Shape); 

D P=zeros(Nsl .NR); 

Mred=s(Index,Index); 

Kred=inv(Mred)*(K 1 l-K12*K22in*K12’); 

Bred=inv(Mred)*(b 1 -K 1 2*K22in*b2); 

Cred_P=Cl_P-C2_P*K22in*K12'; 

Dred_P=D_P+C2_P*K22in*b2; 

[Sae,Ar]^eig(Kred); 

Br=Sae'*Bred; 

CrP=CredP * Sae ; 

Dr_P=Dred_P; 

dispC Default damping value is 1 %'); 

SSS=input( f Do you want to change the damping value for each mode (l=yes,0=no) '); 
Damp=2*sqrt( Ar)*eye(NMnew,NMnew)* .0 1 ; % Percent Damping 

% 

% Modifications to add damping to individual modes 

% L.G. Horta 6-21-00 

% 

ifSSS =1; 

DampValue=zeros(NMnew, 1 ); 
for III=l:NMnew; 

disp([' Frequency - num2str(sqrt(diag(Ar(IH,III)))/(2*pi)) 1 (Hz)']) 
DampValue(IlI)=input('Enter damping value in % = '); 
end 

Damp=2*sqrt(Ar)*diag(DampValue)/100; 

end 


A=[zeros(NMnew,NMnew) eye(NMnew.NMnew); -Ar -Damp]; 
B=[zeros(NMnew,NR); Br]; 

C_l=[Cr_P zeros(Nsl.NMnew)]; 

C_2=[-Cr_P*Ar -Cr_P*Damp]; 

D_l=Dred P; 

D_2=Cr_P* Br+Dr_P ; 


C = [C_1;CJ]; 

D=[D_1;D_2]; 

% 

% Output Definition 
% 


ModelDesc=strvcat([' Number of output states is equal to ' num2str(Nsl*2)],... 

[' The first ' num2str(Nsl) ’ are translations’],... 

[’ From ’ num2str(Nsl+l) ’ to ’num2str(Nsl*2) ’ are 

accel. ’]); 
disp(ModelDesc) 
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